;-------------------------------------------------
; basic parameters for base state density profile
; must agree with the values taken in the eulag
; version used for the corresponding simulation 
;-------------------------------------------------
; add:
; 
; they are naturally configured to the eulag 
; values by defining them with the initial
; values of strat.dat
; --> still need to reset the value of "g"
;_________________________________________________

rh00=rho_strat(0)
t00=temp_strat(0)
th00=t00
;g=1.
if g eq !Null then begin
 print, 'please set the value of the g aceleration ...'
endif
rg=0.4
cp=2.5*rg
cap=rg/cp
capi=1./cap
ss=g/(cp*th00)
press00=rg*rh00*t00
gam=5./3.

;-------------------------------------------------
; Radial discretization -> must agree with 
; /Downloads/code/rk4PACHO.f of the eulag version
; used for the corresponding simulation.
;-------------------------------------------------

zmin=z_strat(0)
zmax=z_strat(-1)
dz=(zmax-zmin)/l
z=indgen(l)*dz+zmin

;-------------------------------------------------
;     ADIABATIC PROFILES !!! 
;
; --> Base state density, pressure, temperature
;     & and theta profiles.
;
; --> Local adiabatic speed of sound profile
;
; --> notice that in this case the theta profile
;     reduces to the case of constant entropy:
;
;     theta_ad = t_ad*((t00*rh00)/(t_ad*rho_ad))^cap
;              = t00    
;
;-------------------------------------------------

rho_ad = rh00*(1-ss*(z))^(capi-1)
p_ad = press00*(1-ss*z)^(capi) 
t_ad = t00*(1-ss*z) 
;theta_ad = t_ad*((t00*rh00)/(t_ad*rho_ad))^cap
theta_ad = t00
cs_ad=sqrt(gam*p_ad/rho_ad)

;-------------------------------------------------
; Compute averages
; 
; To compute the averages in time we take the 
; last "nslices" slices up to the slice "topslice"
; for all the averages  
;-------------------------------------------------

print, 'averages are calculated within the range [topslice,topslice-nslices]'
if nslices eq !Null then begin
 print, 'ERROR ...'
 print, 'please set a value for nslices '
endif else begin
 print, 'current value of nslices = ', nslices
 if topslice eq !NULL then begin
  print, 'ERRROR ...'
  print, 'please choose a topslice to compute the averages'
 endif else begin
  print, 'current value of topslice = ', topslice
  ;________________________________________________________________
  ; SPATIAL AVERAGES -> gives variables of the type f=f(z) only
  ;________________________________________________________________ 
  ;
  ; the total temperature is calculated following the aneslastic
  ; approximation: 
  ;                  (T/T_ad) ~ (theta/theta_ad)
  ; 
  ; where T_ad and theta_ad are the adiabatic profiles.
  ; 
  ; Here we consider two types of temperature fluctuations DT_ad
  ; and DT_e, which are respectively:
  ;
  ;                    DT_ad = T_ad - T      (T_ad -> adiabatic)
  ;                    DT_e  = T_e  - T      (T_e  -> strat.dat)
  ;
  ; The total temperature fluctuations  are correspondingly:
  ;
  ;                    DT_ad/T_ad = (T_ad - T)/(T_ad)
  ;                    DT_e/T_e   = (T_e  - T)/(T_e)
  ;________________________________________________________________
  
  ; "y" averages:
  ; ... for the terms DT/T ....
  
  theaver=total(the,1)/m                        ; <theta>

  ttot_aver=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  for ii=0,n_elements(theaver(0,*))-1 do begin
   for jj=0,n_elements(theaver(*,0))-1 do begin
    ttot_aver(jj,ii)=(t_ad(jj)/theta_ad)*theaver(jj,ii)     
                                                ; <T>=(T_ad/theta_ad)*<theta>
   endfor
  endfor 

  dt_ad=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  dt_e=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  for ii=0,n_elements(theaver(0,*))-1 do begin
   for jj=0,n_elements(theaver(*,0))-1 do begin
    dt_ad(jj,ii)=t_ad(jj)-ttot_aver(jj,ii)      ; DT_ad   
    dt_e(jj,ii)=temp_strat(jj)-ttot_aver(jj,ii) ; DT_e
   endfor
  endfor

  ; "y" averages: 
  ; ... for the term <w2>/cs^2 ....
  
  w2_aver=total(w2,1)/m

  ; "y" averages:
  ; --------- >      Fconv_ad = cp*rho*<DT_ad*w>
  ; ... for the term Fconv_ad/(rho*cs^3) ....
  ; --------- >      Fconv_e = cp*rho*<DT_e*w>
  ; ... for the term Fconv_e/(rho*cs^3) ....
 
  w_aver=total(w,1)/m

  dtw_ad=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  dtw_e=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  for ii=0,n_elements(theaver(0,*))-1 do begin
   for jj=0,n_elements(theaver(*,0))-1 do begin
    dtw_ad(jj,ii)=dt_ad(jj,ii)*w_aver(jj,ii)           ; DT_ad*w
    dtw_e(jj,ii)=dt_e(jj,ii)*w_aver(jj,ii)             ; DT_e*w
   endfor
  endfor

  ;_______________________________________________________________
  
  thew_aver=total(thew,1)/m                            ; <theta*w>
  tpw_aver=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))  
  for ii=0,n_elements(theaver(0,*))-1 do begin
   for jj=0,n_elements(theaver(*,0))-1 do begin
    tpw_aver(jj,ii)=thew_aver(jj,ii)*(t_ad(jj)/theta_ad); <Tp*w>
   endfor
  endfor

  ;_______________________________________________________________

  ; "y" averages:
  ; ... for the term Fkin = 0.5*rho*<u2*w> ....  

  vel2=u2+v2+w2
  tu2w=vel2*w
  tu2w_aver=total(tu2w,1)/m                            ; <u2*w>
 
  ;________________________________________________________________
  ; TIME AVERAGES -> gives a time average over "nslices" time steps
  ; of the previous spatialy averaged variables  of the type f=f(z) 
  ;________________________________________________________________
  ; 
  ;   <<f(z)>>=SUM_[topslice-nslices:topslice] {f(z)} /nslices  ;
  ;________________________________________________________________

  tsize=n_elements(theaver[1,*])            ; last slice
  tmin=(topslice-1)-nslices                 ; first slice for the
                                            ; time average
  
  ;________________________________________________________________
  ; <<Tp*w>>
  tpw_prom=total(tpw_aver[*,tmin:topslice-1],2)/double(nslices)
  ;________________________________________________________________
  
  ; <<DT_ad*w>>
  dtw_ad_prom=total(dtw_ad[*,tmin:topslice-1],2)/double(nslices)
  ; <<DT_e*w>>
  dtw_e_prom=total(dtw_e[*,tmin:topslice-1],2)/double(nslices)
  ; <<u^2*w>>
  tu2w_prom=total(tu2w_aver[*,tmin:topslice-1],2)/double(nslices)
  ; <<DT_ad/T_ad>>
  delta_ad=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  for ii=0,n_elements(theaver(0,*))-1 do begin
   for jj=0,n_elements(theaver(*,0))-1 do begin
    delta_ad(jj,ii)=dt_ad(jj,ii)/t_ad(jj)
   endfor
  endfor
  delta_ad_prom=total(delta_ad[*,tmin:topslice-1],2)/double(nslices)
  ; <<DT_e/T_e>>
  delta_e=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  for ii=0,n_elements(theaver(0,*))-1 do begin
   for jj=0,n_elements(theaver(*,0))-1 do begin
    delta_e(jj,ii)=dt_e(jj,ii)/temp_strat(jj)
   endfor
  endfor
  delta_e_prom=total(delta_e[*,tmin:topslice-1],2)/double(nslices)
  ; <<w^2>>
  w2_prom=total(w2_aver[*,tmin:topslice-1],2)/double(nslices)

;-------------------------------------------------
; Compute  fluxes
;-------------------------------------------------

  fconv=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  fconv_ad=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*))) 
  fconv_e=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  fkin=make_array(n_elements(theaver(*,0)),n_elements(theaver(0,*)))
  for ii=0,n_elements(theaver(0,*))-1 do begin
   for jj=0,n_elements(theaver(*,0))-1 do begin
    fconv(jj,ii)=cp*rho_ad(jj)*tpw_aver(jj,ii)
    ;fconv_ad(jj,ii)=cp*rho_ad(jj)*dtw_ad(jj,ii)
    ;fconv_e(jj,ii)=cp*rho_ad(jj)*dtw_e(jj,ii)
    fkin(jj,ii)=0.5*rho_ad(jj)*tu2w_aver(jj,ii)
   endfor
  endfor

;-------------------------------------------------
; Compute TIME_AVERAGED fluxes
;-------------------------------------------------

  fconv_prom=cp*rho_ad*tpw_prom
  ;convective flux (adiabatic referenced T fluctuations)
  ;fconv_ad_prom=cp*rho_ad*dtw_ad_prom
  ;convective flux (polytropic referenced T fluctuations)
  ;fconv_e_prom=cp*rho_ad*dtw_e_prom
  ;kinetic flux
  fkin_prom=0.5*rho_ad*tu2w_prom

 endelse
endelse
end
